Josephson current through a nanoscale magnetic quantum dot 
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We present theoretical results for the equilibrium Josephson current through an Anderson dot 
tuned into the magnetic regime, using Hirsch-Fye Monte Carlo simulations covering the complete 
crossover from Kondo-dominated physics to n junction behavior in a numerically exact way. Within 
the 'magnetic' regime, l//r > 1 and ea/Y < 1, the Josephson current is found to depend only on 
A/Tk, where A is the BCS gap and Tk the Kondo temperature. The junction behavior can be 
classified into four different quantum phases. We describe these behaviors, specify the associated 
three transition points, and identify a local minimum in the critical current of the junction as a 
function of A/Tk- 



PACS numbers: 74.50.+r, 72.15.Qm, 75.20.Hr 

Recent advances in nanoscale manipulation and fab- 
rication call for a deeper understanding of the effect of 
electronic correlations. Due to the complexity of its the- 
oretical treatment, the interplay between superconduc- 
tivity and magnetism belongs to the least understood 
phenomena in that respect. Here we study the Joseph- 
son current I((f>) through a correlated nanoscale quantum 
dot contacted by s-wave BCS superconductors. At low 
enough temperatures, such a dot is generally described 
by the Anderson impurity model indicated in Fig.^ We 
consider the regime U/T ^> 1 and eo/T -C — 1 , where the 
dot effectively has single occupancy and thus represents a 
spin- 1/2 degree of freedom. Then a complicated interplay 
between this magnetic impurity and the superconductiv- 
ity in the leads sets in. Some aspects of this physics were 
recently observed in Andreev conductance measurements 
for a short multi-wall nanotube 0, 0] • A similar setup 
should also allow to probe the Josephson current in the 
near future, where the ratio A/Tk is widely tunable via 
a backgate voltage. 

In this paper, we provide a detailed analysis and classi- 
fication of all possible phases expected in such an exper- 
iment. We find that only one 'master' parameter A/Tk 
governs this problem, where 



T, 



1 



K 



/ f(7exp[7reo(eo + C/)/r(7] 



(1) 



is the Kondo temperature for normal leads For 
A/Tk <C 1, the Kondo effect survives and is only weakly 
affected by superconductivity, while for A/Tk ^ 1, per- 
turbation theory in T yields an inverted Josephson rela- 
tion I{<j>) = —I c sin0 |1,|1,SI3; where </> = tt represents a 
minimum of the junction free energy F((f>). Such a 7r junc- 
tion behavior was recently reported in Nb-Cu^Nii^-Nb 
systems Q, is related to subgap (Andreev) bound states 
0, Il0j| , and implies broken time-reversal symmetry. In 
both limits, analytical expressions [|| are reproduced by 
our method below. For a magnetic impurity, the Joseph- 
son relation is generally replaced by a more complicated 
dependence on 0. A classification into four types of junc- 
tions, labeled as 0, 0', tt' and tt, follows from the respec- 
tive stability of the (f> — and <fr = 7T configurations 0. 



For a (tt) junction, only <j> = ((f) = it) is a minimum 
of F(4>). For the two other cases, both <fi = 0, tt are lo- 
cal minima, and depending on whether (f> = ((f) = tt) 
is the global minimum, one has a 0' (tt') junction. Us- 
ing I((f>) = (2e/H)dF((p)/d(f), the phase boundaries can 
be directly read off from the I((f>) curves. For instance, 
the 0'-tt' transition point is determined by the condition 

= 0. 

The theoretical description of the resulting phase dia- 
gram is difficult, and despite intense efforts over the past 
few years 0, E3, E3, E3, E3, E3, ES EH E3, E^l , a satisfac- 
tory physical picture has not been obtained so far. This 
paper provides numerically exact results for a magnetic 
Josephson junction from Monte Carlo (MC) simulations, 
using a Hirsch-Fye algorithm 0, [2(], El adapted to this 
problem. We find that previous approximate theories 
for this problem, relying on the non-crossing approxi- 
mation (NCA) for U -> oo [H El El, mean-field ap- 
proaches 0, , perturbative schemes |l4 EH lliil. Il7l | , or 
the numerical renormalization group (NRG) |l2tll8j ]. lead 
to incomplete and sometimes even qualitatively inaccu- 
rate predictions. Although the existence of the above- 
mentioned phases follows already from mean-field theory 
, their respective stability and the actual phase bound- 
aries have not been reliably established. Moreover, the 
critical current defined by 



I C (A/T K ) = max \I(</>,A/T K ) 



(2) 



has a non-monotonic behavior as a function of A/Tk 
which has been missed by all previous studies. The pre- 
dicted minimum as well as the phase boundaries specified 
in Eqs. (|8I9I1U|I below should be observable in state-of- 
the-art experiments. 

We stud y th e 'canonical' model for this problem, see 
Refs. [1 1 M ElElElElElElElEland Fig. □ 
For a symmetric situation, 
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Nambu indices [i,v — 1,2 and time indices and 



BCS 

A, (f>/2 



FIG. 1: Schematic setup of an Anderson dot with local energy 
eo and charging energy U > coupled to BCS leads. For 
simplicity, we assume identical BCS gap A on both sides, 
with phase difference <j> across the dot and hybridization V. 



H = 



k,p a 

(Ae^/ 2 c{ !p (fc)c{ jp (-fc)+h.c. 



(3) 



[ £ o + — J K + H) ~ 2"( n T - n i) > 



where c ffjJ ,(fc) is the electron operator for lead p = L/R = 
± and spin a =|, |, with single-particle dispersion e k - 
Moreover, n CT = with the dot electron operator d a , 
and T = 27rp |^| 2 for lead density of states p and hop- 
ping matrix element P. In all simulations reported below, 
T/A = 0.1 (we put h = ks = 1), which corresponds to 
a temperature T ~ 100 mK for the setup of Refs. PIQ- 
By comparing to analytical results for A/Tk 3> 1 and 
<C 1, this appears to be quite close to the ground-state 
limit. 

To formulate the MC scheme, we construct the 
imaginary-time path integral representation under this 
model for the Joscphson current. Discretizing imaginary 
time in steps of size 6=1 jPT for P discretization points, 
the last term in Eq. © can be decoupled by auxiliary 
Ising spins s& = ±1 defined at times r k = kS, where 
k = 1,...,P, using the discrete Hubbard-Stratonovich 
transformation 191 



3 (C/5/2)(n T -n i ) 2 (r fc ) 



1 E ■ 

2 <H 

s*=±l 



-As fc (n T -n i )(r fe ) 



(4) 



where A = cosh" 1 [cxp({7 8/2)]. All fermions characteriz- 
ing the dot and the leads are then free and can be inte- 
grated out. Thereby the Josephson current is expressed 
in terms of a cyclic ID Ising spin chain with non-standard 
long-ranged spin-spin interactions, 



J2 M det(G" x + As) Tr j [G" 1 + As] _1 E J | 



E W det(G-i+As) 



(5) 



where 1$ = eA/h is the critical current in the uni- 
tary limit, s has matrix elements s-y.iw = SiS^S^ with 



Gr. 1 = 5 2 TV e -**>n*(i-J) 

r 



-iuJ n T + (e + U/2)t 3 



The Joscphson self energy is 



(iuj n T + Acos((/)/2)t 1 ) 



(0) 



p — iuj n S(i— j) 

t J = -( ( 5T) 2 rsin(0/2)r 1 V - (7) 

All summations over fcrmion Matsubara frequencies 
Lu n = (2n+ 1)ttT are restricted to \uj n \ < tt/S, since finite 
S implies the existence of a UV cutoff. The Pauli matrices 
r, act in Nambu space, with tq = diag(l, 1). We mention 
that a modified version of this approach could also access 
the case of unconventional superconductor leads, where 
additional features are expected [22L l23l| . 

Equation JSJ allows to compute the equilibrium 
Josephson current for arbitrary parameters under a MC 
scheme. Finite-<5 corrections can be eliminated by ex- 
ploiting the theorem that such corrections must scale 
cx S 2 for any Hermitian observable pp| |. For given phys- 
ical parameters, we thus compute I(<f) for a few (small) 
values of 8, and then perform the extrapolation 8 — > us- 
ing a linear-regression fit. The S 2 scaling is well obeyed 
for AS < 0.1 even for U/A = 20, T/A = 5, leading to 
discretization numbers P « 140 to 270. Then no system- 
atic errors are present, and numerical results are exact 
within stochastic MC error bars. Due to the absence of 
particle-hole symmetry, this algorithm has a sign problem 
for A =^ 0, which manifests itself in occasionally negative 
determinants in Eq. (JSJ) ■ Fortunately, this problem is very 
weak in our implementation, with average sign above 0.7, 
and hence does not restrict the method in practice. Lo- 
cal flip updates of the Ising spins under the standard 
Metropolis algorithm were sufficient to ensure rapid equi- 
libration and satisfactory MC acceptance rates. Under a 
spin flip s k — > — Sjt, corresponding to As — > Xs + w, where 
^■p = —2XskSikSjkS^, the change in weight is given 
by the ratio R of new and old determinants appearing in 
Eq. ©, which can be evaluated analytically. With the 
matrix D = (G _1 + As) -1 , we find 

R = det(l + Dw) = 1 - 2Xs k (D kk>n + D kU . 22 ) 

+AX 2 [Dkk ,11-Dfefe ,22 — ffefe,12-Dfefe,2l]- 

Thereby the costly explicit calculation of determinants is 
avoided. Similarly, D can be updated without explicit 
matrix inversion. Typically, 10 6 MC samples were accu- 
mulated to obtain each data point below. On a 2 GHz 
Xeon processor, our code performs at a speed of 3.7 CPU 
hours per 10 5 samples for P = 180. 

First, the code was checked against analytical solutions 
available for small and large A/Tk which were accu- 
rately reproduced, see, for instance, the inset of Fig.0for 
A/Tk 1- Our simulations are therefore able to cover 
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the complete crossover from Kondo-dominated physics 
to a 7r junction. (The case A = has also been studied 
in Rcfs. [Uliil 

using this method.) Moreover, we have 
checked for many parameter sets that universality is ful- 
filled, i.e., taking different values for eo/A, U/A and/or 
r/A only affects the Josephson current via the corre- 
sponding change in A/Tk- However, the two curves for 
€q = —U/2 shown in the main part of Fig.[3]dcmonstrate 
that profound differences can arise even for the same 
A/TV, if U/T is not sufficiently large. We found U/T > 5 
necessary to ensure universality; otherwise charge fluctu- 
ations may alter the Josephson current even qualitatively, 
see Fig. [5] Universality then holds only in the true mag- 
netic regime, and not away from it. We mention in pass- 
ing that the critical current never exceeded Iq, in contrast 
to the prediction in Ref. [Tlj . 

We then continue by discussing the full crossover be- 
tween these two limiting cases. Numerical results shown 
below were obtained for 5 < U/T < 10 and 0.1 < A/T < 
10, putting e = —U/2. In Figures and H| repre- 
sentative data for I(4>) covering this crossover are pre- 
sented. Furthermore, in Fig. [3] we show the critical cur- 
rent J2J) as a function of A/Tk- For very small A/Tk, 
the Kondo effect is dominant, and we find a junction 
with a non-sinusoidal current-phase relation close to the 
unitary limit, I c = 2o 6|. For larger A/Tk, the rela- 
tion is more sinusoidal again, and with increasing A/Tk, 
the critical current I c decreases. Moreover, above a first 
transition point 



(A/TkV =2.8 ±0.1, 



(8) 



for 4> slightly below n, the current is negative, and un- 
der the above classification scheme, we thus enter the 
0' phase, see Fig. O By increasing A/Tk further, one 
eventually reaches a second transition point at 



(A/7VV 



7.2 ±0.2. 



(9) 



At the transition point, f d<j) I(4>) = 0. We have now 
reached the tz' phase, where <j> = ir is already the global 
minimum, but = still represents a local minimum. 
The true 7r phase, with <f> = tt as the only minimum, is 
eventually reached at 



(A/Tk) 11.0 ±0.3. 



(10) 



For A/Tk 3> {A/Tk)tt'tt, the inverted Josephson relation 
/ = —I c sin <f> valid in the deep tt junction limit [6j is 
finally recovered, sec Figs. [21 and 

The transition points reported here are at quite dif- 
ferent locations than thought previously. For instance, 
NRG calculations 0, find that the transition into 
the 7r phase occurs at A/Tk ~ 2.4, where according 
to our data the junction is still in the phase. The 
differences to NCA and/or mean- field results are even 
more drastic. According to our simulation data, the 7r 
phase covers a much smaller region in parameter space 
than thought previously, while the intermediate 0' and 
7r' phases extend over a significant range in A/Tk, 
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FIG. 2: MC results for the Josephson current in units of 7o = 
eA/h for two parameter sets with A/Tk ~ 23. Squares are 
for U/A = 20, r/A = 3.44, and circles for U/A = 4, T/A = 1. 
Unless noted otherwise, in all figures, curves are guides to the 
eye only. Stochastic MC errors are always smaller than the 
symbol size or indicated by vertical bars. Inset: Current- 
phase relation in the deep tt junction regime. MC results for 
A/Tk = 6.5 x 10 4 coincide with the analytical result (dashed 



o A/T K =3.0 
□ A/T K =4.6 
o A/T K =6.9 




FIG. 3: MC data for the current-phase relation at small-to- 
intermediate A/Tk- 



and therefore should be readily observed in practice. 
Remarkably, the critical current shown in Fig. [S] displays 
non-monotonic behavior. From the analytically known 
limits, I C {A/T K -> 0) -> I and I C (A/T K -> oo) -► 0, a 
naive guess is to expect that I c (A/Tk) just drops moiio- 
tonically with increasing A/Tk- Such a behavior is in 
fact predicted by previous work, see, e.g., Ref. n^j. Our 
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o A/T K = 6.9 
A/T K = 9.7 
A/T K = 11.4 
A/T K =379 




FIG. 4: Same as Fig. [3] for intermediate-to-large A/Tk- 
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simulations point to a more complicated picture, where 
I c (A/Tk) is characterized by a local minimum. This 
minimum occurs at 



8.2 ±0.2, 



[11) 



which is close to but above the transition point JJjJ for the 
O'-tt' transition. After reaching a local maximum slightly 
below the tc'-tt transition, the critical current then drops 
monotonically throughout the 7r regime. At first sight, 
the local minimum in I c appears to resemble the exper- 
imentally observed oscillations in the critical current as 
a function of cither temperature or length of the junc- 
tion ||. However, while those oscillations can be traced 
back to Andreev bound state crossings, such a simple, es- 
sentially mean-field type reasoning does not apply here, 
see also Ref. for a closely related discussion. The ap- 
pearance of a local minimum in I c thus indicates a subtle 
many-body effect. 
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To conclude, we have presented a numerically exact 
study of the Josephson current through a nanoscale dot. 

° In the magnetic regime, U/T 3> 1 and eq/T < 1, our re- 
sults reveal a rather complex behavior that is governed by 

1000 A/Tjc as the only tuning parameter. These predictions 
should be observable in experiments on short carbon nan- 
otube quantum dots. 



FIG. 5: MC results for the critical current Q as a function 
of A/Tk- The arrows mark the boundaries between different 
phases (see text). The inset shows the region around the 
minimum on a magnified scale. Error bars are always smaller 
than symbol size. 
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